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Abstract. We present a re-reduction of archival CCD frames of the doubly imaged quasar 0957+561 using a new 
photometry code. Aperture photometry with corrections for both cross contamination between the quasar images 
and galaxy contamination is performed on about 2650 _R-band images from a five year period (1992-1997). From 
the brightness data a time delay of 424.9 ± 1.2 days is derived using two different statistical techniques. The 
amount of gravitational microlensing in the quasar light curves is briefly investigated, and we find unambiguous 
evidence of both long term and short term microlensing. We also note the unusual circumstance regarding time 
delay estimates for this gravitational lens. Estimates by different observers from different data sets or even with the 
same data sets give lag estimates differing by typically 8 days, and error bars of only a day or two. This probably 
indicates several complexities where the result of each estimate depends upon the details of the calculation. 
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1. Introduction 

The first reported example of gravitational lensing, the 
twin quasar QSO 0957+561, was discovered in 1979 by 
Walsh et al. (1979). It is one of the most studied objects in 
modern cosmology, and the research and monitoring cam- 
paigns have mainly been fueled by the desire to measure 
the time delay, and thereby, to get an independent and di- 
rect estimate of the Hubble parameter (Refsdal 1964). In 
addition, several groups have tried to analyze the extrinsic 
variability in the light curves. This variability is assumed 
to be caused by gravitational microlensing (ML) by stars 
or MACHOs in the lensing galaxy (as predicted by Chang 
& Refsdal 1979). 

QSO 0957+561 is a doubly imaged quasar at a redshift 
z = 1.41, with the components A and B separated by 
about 6'.'2. 1 A massive, elliptical cD galaxy (named Gl) 
at z = 0.36, located only ~ 1" from the center of the B 
image, seems to be the principal lensing object. 

The closely juxtaposed quasar images and the ex- 
tended brightness profile of the lens galaxy make accu- 
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1 The literature consistently quotes 6'.'l, although accurate 
HST and VLBI astrometry yields values of 6'.'169 and 6'.'175, 
respectively (Bernstein et al. 1997). 



rate photometry a challenge. During the 1980s and mid- 
1990s, standard aperture photometry was performed with- 
out any corrections for the light contamination between 
the quasar components (crosstalk) or from the lens galaxy, 
see e.g. Schild & Cholfin (1986); Schild (1990); Kundic 
et al. (1995). Later reduction schemes have tried to ad- 
dress the above-mentioned problems in order to reduce the 
chance of correlated (seeing-dependent) brightness varia- 
tions in the light curves; e.g. Colley & Schild (1999, 2000); 
Serra-Ricart et al. (1999). Precise photometry is neces- 
sary for time delay determinations and for investigations 
of possible microlens-induced fluctuations in the bright- 
ness records. 

In spite of extensive observations by several groups, 
the time delay (r) between the two quasar images has 
proved hard to determine. Complicating factors include 
heterogeneous data sets, large temporal gaps in the data 
sets, and additional variability in one or both of the quasar 
images (microlensing). Even more than 15 years after the 
discovery, the time delay was not determined. However, 
there were two favored candidates; ~ 540 days and ~ 415 
days. The results of the different investigations prior to 
1997 are summarized in Haarsma et al. (1997). Kundic 
et al. (1997) convincingly settled the long-standing con- 
troversy in favor of the lower value, finding r = 417+3 
days. Since 1995 different groups have reported values of r 
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in the range 416-425 days. The results seem again to con- 
centrate around two values; 417 days (Kundic et al. 1997; 
Pelt et al. 1998; Colley & Schild 2000) and 424 days (Pelt 
et al. 1996; Oscoz et al. 1997; Pijpers 1997; Serra-Ricart 
et al. 1999; Oscoz et al. 2001). 

QSO 0957+561 was the first system to provide strong 
indications of microlensing effects; uncorrelated bright- 
ness variations between the A and B images were found 
by Vanderriest et al. (1989). Several researchers have re- 
ported microlens-induced variability in the quasar light 
curves. Pelt et al. (1998) found unambiguous evidence 
of long time scale (order of several years) microlensing 
in the "difference light curve" (DLC; A light curve mi- 
nus time-shifted B curve). Results are ambiguous when 
it comes to the short time scale (lasting a few months) 
and rapid (less than a few weeks) microlensing events. 
Schild & Thomson (1995a), Schild (1996) and Colley & 
Schild (2000) have reported interesting high-frequency fea- 
tures in the brightness record, having amplitudes of only 
0.03 — 0.05 mag and time scales of months and even weeks. 
Goicoechea et al. (1998) also found fluctuations which 
could be associated with microlensing events. However, 
Schmidt & Wambsganss (1998), Wambsganss et al. (2000) 
and Gil-Merino et al. (2001) all found DLCs with no clear 
microlensing signature, and notably no short time scale 
events with |Am| > 0.05 mag were observed. Gil-Merino 
et al. actually showed that the fluctuations in their DLC 
could be due to (several) observational noise processes. To 
reveal any rapid fluctuations caused by microlensing, high 
quality images with good temporal sampling are required. 

This paper is mainly a summary of some results from 
a Master's thesis project by Ovaldsen (2002) undertaken 
at the Institute of Theoretical Astrophysics, University 
of Oslo, Norway. We shall here concentrate on the aper- 
ture photometry scheme, the time delay estimation and 
microlensing investigation. The data set consists of some 
2650 archival CCD images of QSO 0957+561 covering a 
period of nearly five years (June 1992 - April 1997). This 
data set has previously been reduced by one of the au- 
thors (RES), but with cruder corrections for crosstalk and 
galaxy contamination. In the next two sections we discuss 
the data set and briefly present the main principles of 
our photometry scheme. Then, from the final A and B 
light curves (Sect. 4), the time delay is determined us- 
ing two different statistical techniques (Sect. 5). In Sect. 6 
we briefly investigate the microlensing residual. The re- 
sults are summarized and discussed in Sect. 7. Our new 
photometry gives, among other things, a time delay that 
differs significantly from the result we obtain when em- 
ploying the same method on the old RES brightness data. 

2. Data set 

RES and collaborators have monitored this lens system for 
over a decade and amassed a large data set. Here we use 
a subset consisting of around 2650 i?-band CCD images 
taken with the 1.2m telescope at Fred Lawrence Whipple 
Observatory atop Mt. Hopkins, Arizona, during a five year 



period from June 1992 to April 1997. About 200 images 
were discarded at an early stage due to various CCD de- 
fects, cosmic ray hits, guiding errors, bad pre-processing 
etc. There are usually 4 or 6 frames per night. 

Although taken with the same telescope, the quality of 
the frames varies considerably Cosmic rays, and especially 
bad pixels and bad columns, occur frequently Quite a few 
frames exhibit varying background levels, not only in the 
form of a gradient across the image, but as bright or dark 
"patches" at certain locations. Such frames may not have 
been properly calibrated. (Other artifacts from poor pre- 
processing are also seen). 

The image headers do not contain all the desired in- 
formation, e.g. the pixel size and the gain factor are often 
missing. The gain is fixed to 2.3 e~/ADU. The pixel size 
is computed empirically for each frame, using the calcu- 
lated positions of typically 6 field stars and the astrome- 
try presented in the Guide Star Catalog II (GSC-II), see 
Table l. 2 

Two different CCDs are employed. The scale of the 
first one is approximately 0.65"/pixel (binned mode) and 
0.32"/pixel (unbinned), and that of the second one is 
0.70"/pixel. The range of seeing values (FWHM) is ap- 
proximately 1" — 5", with a mean value of around 2". The 
global background is mainly between 100 and 2000 ADU 
(92% of the frames). The stellar images are typically non- 
circular, the PSF having a mean cllipticity of 0.09, equiva- 
lent to an axis ratio of 1.1. We also note that the PSF often 
departs from elliptical symmetry. The coma-like appear- 
ance is probably due to tracking errors and astigmatism 
in the camera optics. 

The sampling of the observations must be regarded as 
very good. Besides the gaps in the summer months, the 
one day interval dominates. More than 90% of all time 
separations between consecutive observation runs are less 
than eight days. 

A very different and more homogeneous data set, com- 
prising some 1000 R- and F-band frames obtained over 
four consecutive nights, is discussed in a forthcoming pa- 
per (Ovaldsen et al. 2003). 

3. Aperture photometry scheme 

The software used to reduce and analyze the CCD frames 
was developed by JT and JEO. The entire package 
was written almost from scratch in the Interactive Data 
Language (IDL) 3 ; it is specially adopted to the 0957+561 
twin quasar system. Several sub-routines had been imple- 
mented by JT when working on other quasar lens sys- 
tems. Many of these remained unchanged. All steps are 
automated; from detection and localization of objects, via 
field star photometry and calibration, to the actual quasar 
photometry We will only describe the main features of 
our photometry scheme. The automatic source detection 

2 The GSC-II is a joint project of the Space Telescope Science 
Institute and the Osservatorio Astronomico di Torino. 

3 A product of Research Systems, Inc. 
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program, background determination, centering algorithms 
etc. are not discussed here. We refer to Ovaldsen (2002) for 
a complete and detailed treatment of the entire package. 

3.1. Field star photometry and calibration 

The separation between the two quasar images of <~ 6" 
motivates the use of 3" radius apertures. Using the same 
aperture size for the comparison stars as for the target 
objects makes it easier to transform the quasar intensities 
into standard magnitudes. We use the stars F, G, H, E, D, 
X and R for the stellar photometry, see Fig. 1 and Table 1. 
Of course, the use of such small apertures is only valid 




Fig. 1. The field surrounding QSO 0957+561A,B. North is up, 
east is left. The seven field stars (F, G, H, E, D, X and R) 
as well as the quasar components (A and B) are indicated. 
The lens galaxy, not visible here, is located approximately 1" 
from B, slightly to the left of the vector from B to A. The 
frame is 4.7 arcmin on each side. (This particular image is 
taken with the 2.4m Hiltner telescope, Michigan-Dartmouth- 
MIT observatory, and is of better quality than the frames in 
our data set both in terms of seeing and spatial resolution; 
pixel size=0.275"/prxel, seeing w 1"). 

if they collect the same fraction of the total light for all 
point sources, equivalent to the assumption that all point 
sources on a frame have the same PSF. We assume that 
this is approximately the case. Preferably, one should ob- 
serve reference stars with the same spectral distribution 
as the primary targets; this would reduce the error when 
calculating the magnitudes of the target objects. In the 
case of 0957+561, the two quasar images are bluer than 
the field stars. However, since we only have single band 
observations, we are not able to correct for any color ef- 
fects. 



Table 1. Reference star astrometry from GSC-II. Both right 
ascension, a, and declination, S, are given in degrees. The right- 
most column quotes the names of the stars as they appear in 



GSC-II. 


Star 


a 


8 


GSC-II id 


F 


150.41117868 


55.90692350 


N212232175 


G 


150.40026155 


55.89505952 


N212232178 


H 


150.36551899 


55.89125437 


N212232182 


E 


150.37040443 


55.87405306 


N212232187 


D 


150.36712201 


55.86871933 


N212232190 


X 


150.41965483 


55.87058001 


N212232189 


R 


150.40374136 


55.86788639 


N21223213199 



The large pixel size (mostly ~ 0.65"/pixel) combined 
with the relatively small apertures (radius=3") give a 
quite irregular polygon on the pixel array. To simulate 
a "perfect" circular aperture we apply a weighting scheme 
for the pixels which lie on the border. The value of a partial 
pixel is calculated as the original pixel value multiplied by 
the ratio of the partial pixel area to the total (square) pixel 
area. We also tried to quantify the implication of a non- 
zero brightness gradient across the aperture border. This 
second-order correction, however, proved insignificant. 

The local background level is calculated from 20 small 
apertures arranged in a circle of radius 20" around the 
object of interest. The apertures containing cosmic rays, 
bad columns, sources etc. are automatically discarded. 

We use seven comparison stars (Fig. 1) to determine 
the calibration level needed to put the quasar magnitudes 
on the standard system. The instrumental intensities are 
compared to the reference values, and any (5a) outliers 
are registered. Some frames contain fewer than seven com- 
parison stars, but we require a minimum of three to pro- 
ceed. If there is more than one outlier, the frame is simply 
discarded. Our fundamental assumption is thus that the 
intensities of all the present stars (except one possible out- 
lier) should be consistent with the reference values. The 
measurement errors are taken into account. With this pro- 
cedure we make sure that the calibration constant is cal- 
culated from "well-behaved" stars and, consequently, that 
the quasar magnitudes will be as accurate as possible. 

All measurement uncertainties (i.e. the standard devi- 
ation of the aperture intensity I) are calculated using the 
formula 



a(I) = ^I/g + n a (l + ^jv b . (1) 

I is the background-subtracted source intensity (in ADU), 
g is the gain factor of the CCD, Vb is the variance of the 
background, and n a and rib is the number of pixels used in 
the determination of the source intensity and background 
level, respectively. The Vb parameter not only measures 
the variance of the sky level, but also fluctuations (in- 
homogeneities) due to faint background sources and the 
CCD readout noise. Hence, Vb is always larger than the 
Poisson variance of the sky level. 
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3.2. Photometry of the two quasar components 

As previously mentioned, we use 3" radius apertures cen- 
tered on each quasar image. These apertures will be sub- 
ject to seeing- (and ellipticity-) dependent light contam- 
ination from the neighboring quasar component and the 
underlying Gl lens galaxy. 

3.2.1. Galaxy subtraction 

To correct for the galaxy's light contamination, we decided 
to subtract from each frame a synthetic model of Gl. Upon 
request, G. Bernstein kindly provided surface brightness 
data obtained from HST observations in the U-band and 
Kitt Peak observations in the i?-band (see Bernstein et al. 
1997). In order to find the color offset for the HST data, 
we simply looked for what gave the best agreement in the 
overlap area of the profiles. The "correction" V — R = 1.3 
seemed to merge the profiles well. Although the study by 
Bernstein et al. (1997) indicated an ellipticity gradient 
and isophote twist, we decided to model the galaxy with 
fixed values for the ellipticity (0.28) and position angle (53 
degrees) - compare with their Fig. 2. 

The position of the center and the orientation of 
the semi-major axis are calculated from the positions of 
the quasar images and from the relative astrometry of 
Bernstein et al. (1997). To synthesize the galaxy we start 
by oversampling the pixels four times. The value of each 
sub-pixel is computed by interpolating the brightness pro- 
file quadratically. The ellipticity and position angle are 
taken into account. When determining the calibration (or 
zero) level for the galaxy, the small 3" apertures do not 
suffice. Gl is an extended object whose profile is much 
broader and totally different from that of the stars. For 
this reason it is important to calibrate a resolved object 
like Gl with the "total" light from the comparison point 
sources, here taken to be the flux in apertures of radius 
12". 

The model is finally "smeared out" in accordance with 
the seeing on each particular image. This is done by 
convolving the synthesized galaxy with the image PSF. 
Having performed the proper scaling and positioning, the 
convolved galaxy image is simply subtracted from the 
frame. 

3.2.2. Crosstalk correction 

Several methods to minimize cross contamination be- 
tween the A and B images were explored, some of which 
were similar to the procedure in Colley & Schild (1999). 
However, because we have to calculate the PSF for each 
frame (used in the modeling of the synthetic lens galaxy) , 
we decided to utilize one of the characterizing features 
of the PSF-fitting technique. The A and B images are 
cleaned from the frame in an iterative fashion, thereby 
allowing aperture photometry to be performed on each 
quasar image after the galaxy is subtracted and after the 
neighboring twin is cleaned from the frame. The cleaning 



works well for a wide range of seeing conditions, and this 
way of eliminating the crosstalk between the A and B im- 
ages proved to be significantly more robust than the other 
methods (it is, for instance, less sensitive to bad columns, 
bad pixels etc.). 

4. Photometric results 

4.1. Stellar photometry 

2486 frames were analyzed with respect to reference star 
photometry, and the calibration procedure (see Sect. 3.1) 
accepted 2028 frames. Fig. 2 shows the magnitudes of the 
seven reference stars for all accepted frames, as a function 
of Julian Day (J.D.). Table 2 shows some statistics for 
each star. We remark that the stars F, G, H, E, D, X and 
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Fig. 2. Aperture photometry of the field stars around QSO 
0957+561. Plotted are the i?-band magnitudes of the reference 
stars from all accepted data frames, as a function of Julian Day. 
Dashed red line is the mean value, dotted green lines are the la 
limit. The stars are arranged in order of decreasing brightness. 
Note that the scaling for the R star is different from the others. 

R were saturated on 280, 289, 97, 0, 2, 719 and frames, 
respectively. 

From each frame we only calculated the magnitudes 
of the comparison stars which were unsaturated and had 
intensities consistent with the reference values. The cali- 
bration constant was computed from the same stars. As 
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Table 2. Statistics of the field star aperture photometry. For 
each star the mean magnitude and standard deviation are tab- 
ulated, along with the mean of the individual formal errors, 
and reference magnitude. 



Star 


Mean _R-mag 


O" mag 


Mean formal error 


Ref. ma 


F 


13.757 


0.0033 


0.0030 


13.757 


G 


13.730 


0.0029 


0.0029 


13.730 


H 


13.977 


0.0040 


0.0032 


13.977 


E 


14.783 


0.0048 


0.0047 


14.782 


D 


14.515 


0.0045 


0.0041 


14.515 


X 


13.422 


0.0032 


0.0027 


13.422 


R 


16.336 


0.0135 


0.0123 


16.335 



we see from the plots, there are a few "outliers" . Some of 
these can be identified as points with large error bars, but 
some are true outliers in the sense that they should have 
been discarded. The calibration constant is not necessarily 
biased by such outliers in all cases, because it is derived 
from several (3-7) stars. 

As expected, the scatter in magnitudes increases for 
fainter stars. The R star is by far the faintest, and conse- 
quently has the largest standard deviation for the calcu- 
lated magnitudes, i.e. 13.5 mmag. Given that its bright- 
ness is ~ 0.3 mag greater than the quasar images, this 
should indicate the minimum general dispersion to be 
expected in the quasar light curves due to photometric 
noise. After all, the photometry of the two quasar images 
is much more complicated than that of an isolated star; 
the galaxy subtraction and the cross talk correction obvi- 
ously increase the error budget of images A and B. 

4.2. Quasar photometry 
4.2.1. Binning and censoring 

The number of data points (or accepted frames) per night, 
n, varies from one to six. The observation were made 
within a small time interval, thus it seems appropriate 
to combine all data points for a particular night, and only 
quote "binned" magnitude values. It is almost impossi- 
ble to perform a rigorous statistical analysis with so few 
data points (and an unknown number of outliers due to 
erratic photometry of the twin images). We decided to 
address this issue in a simple and transparent way. From 
the image A magnitudes and the image B magnitudes on 
each night, we computed the corresponding median val- 
ues. With this approach, single outliers do not bias the 
results too much, at least for n > 3. Obviously, for n = 
1 or 2, the median- filtering will not throw away possible 
outliers. However, without any a priori knowledge, this 
is about the best we can do. We quote error bars which 
correspond to the median points. 

Before the median filter was applied we censored the 
data, accepting only images with seeing < 3" FWHM, 
background level < 3000 ADU, PSF axis ratio < 1.3 and 
AB-separation 6'.' 175 ± 0'.'05. The final data set was then 



reduced to 1720 images. The subsequent "binning" yielded 
422 data points for each quasar image. 

4.2.2. Crosstalk and galaxy contamination 

Aperture photometry on images A and B was performed 
both with and without the corrections for crosstalk and 
galaxy contamination (see Sect. 3.2), so that we could 
check the performance. We now make a few comments 
on the results from analyzing the 1720 accepted images. 

Fig. 3 displays how the light from the two quasar im- 
ages affects each aperture. We emphasize that this effect 

Aperture cross talk 



-1 r 




i.O 1.5 2.0 2.5 3.0 

FWHM seeing (") 

Fig. 3. The percent change in the A and B aperture flux when 
correcting for crosstalk, as a function of seeing. The crosstalk 
from the B (A) image into the A (B) aperture is marked with 
open squares (triangles). The corrections are virtually the same 
for the two apertures. A best fit quadratic curve is overplotted. 

is not only a function of seeing; the ellipticity is certainly 
also a factor. In particular, consider an image where the 
PSF is highly elliptical and has a semi-major axis parallel 
to the line joining the centers of A and B. The crosstalk 
would be larger here compared to the case where the semi- 
major axis is perpendicular to the AB vector. (In fact, the 
scatter of the corrections can be significantly reduced by 
imposing stricter limits on the ellipticity). 

For completeness, a least squares second order poly- 
nomial fit was computed using all the data points (see 
also Colley & Schild 2000). The formula for the intensity 
corrections, SI, reads 

51(a) = [0.944 - 0.546s + 0.397s 2 ] % , (2) 

where s = FWHM ("). As can be seen in the figure, the 
curve fits the data quite well. The scatter is partly due 
to ellipticity, but some of the deviant points are a result 
of bad cleaning of the quasars images; on some frames 
the flux in the measuring aperture is still affected by the 
neighboring quasar image, which has not been properly 
removed/subtracted. 

We also measured the flux in the quasar apertures 
before and after the galaxy model had been subtracted. 
Fig. 4 shows the light contribution from the galaxy to 
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the A and B apertures as a function of seeing. Note that 
crosstalk between A and B has already been "eliminated" . 
As seeing deteriorates, galaxy light systematically seeps 



Galaxy contamination In quasar apertures 
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FWHM seeing (") 
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Fig. 4. Light contamination (in percent) in the A and B aper- 
tures due to the Gl lens galaxy, as a function of seeing. Best 
fit quadratic curves are also plotted. 

out of the B aperture, but into the A aperture. Best fit 
quadratic curves are overplotted to guide the eye. Some of 
the scatter in the plots is due to the fact that A and B itself 
varied during the observational period (the Gl contribu- 
tion is compared to the A and B fluxes on each particular 
frame.) The contribution to the A image aperture is be- 
tween 3 and 5%, and has a moderate scatter. The B image 
aperture has corrections of roughly 20—18% due to the Gl 
galaxy. Here, the scatter is rather large, having a "full am- 
plitude" of ±2%. It does not seem to increase with seeing. 
This probably indicates that the calibration/zero level for 
the galaxy is determined equally well (or poor!) for the 
whole range of seeing conditions. Comparable corrections 
for subtracting the lens galaxy contribution were discussed 
in detail by Colley & Schild (2000), whereas the original 
RES reductions incorporated subtraction of a fixed 18.34 
magnitude correction for the lens galaxy contribution to 
the B aperture flux. 

4.2.3. Quasar light curves 

Fig. 5 displays the light curves of QSO 0957+561A,B cor- 
responding to the period June 1992 to April 1997. 4 We 
note that the light curves show variability on both short 
(order of weeks) and long (order of years) time scales. For 
some periods there is also an apparent "zero lag" correla- 
tion between A and B. This is best seen for J. D. -2448000 
> 1700. Such a correlation should in general not exist, be- 
cause the signal from image B lags A by some r ~ 420 

4 The data table is available at 

http : //www. astro .uio . no/~jeovalds/DQmags .html. 



days. We have not been able to identify the cause of this 
frame-by-frame correlation. It has also been reported by 
Colley et al. (2003), and is always presumed to be some 
systematic effect caused by errors in the photometry; how- 
ever the amplitude exceeds any known error source. 

Although our "binning" scheme only uses the nightly 
median magnitude value for each quasar image, we can 
still estimate the dispersion for nights with two or more ac- 
cepted frames. The mean standard deviations of the mag- 
nitudes on each night are 12 mmag and 11 mmag for A and 
B, respectively. We note that the mean of the formal error 
bars, as seen in Fig. 5, is 17 mmag for both quasar images. 
The formal error bars are rather conservative, as they in- 
clude the formal errors from Poisson statistics, galaxy sub- 
traction and calibration (see Ovaldsen 2002 for details). 

We also made a rough and simple estimate of the day- 
to-day dispersion within the A and B brightness data: 
First each light curve was smoothed with a 7-point me- 
dian filter (making sure not to filter over gaps greater 
than three days). Then the original data (A or B) was 
subtracted from the corresponding median-filtered curve. 
We allowed a maximum time gap of 1.5 days between 
two data points to be subtracted. The residuals should 
thus probe fluctuations in the A and B light curves on 
this time scale. We assume that this very short time scale 
variability is not dominated by microlensing effects. The 
standard deviations of the residuals are cr^ sld = 11 mmag, 
and <7g sld = 10 mmag. These values are quite consistent 
with similar estimates for the image set made by Schild & 
Thomson (1995b), who found 9.5 and 12.0 mmag; however 
their reductions lack the corrections for aperture crosstalk 
and galaxy subtraction that are strictly functions of see- 
ing, and are more susceptible to systematic errors on time 
scales relevant to seeing changes. 



5. Time delay 

The main analysis is performed using a method based on 
dispersion estimates, but we also explore a different tech- 
nique based on \ 2 minimization. We use the data corre- 
sponding to Fig. 5. 

5.1. Dispersion estimates 

The algorithm for the Dispersion estimation technique is 
included in the ISDA (Irregularly Spaced Data Analysis) 
package, designed by J. Pelt to perform various tasks on 
irregularly spaced time series. It is discussed extensively 
in Pelt et al. (1994, 1996), so we provide here only a short 
review. The principle of the Dispersion method is simply 
to measure the dispersion, D 2 , between the A and B image 
light curves for different time shift values, r. The true 
value should show up as a minimum in the dispersion 
spectrum, D 2 (t). By dispersion we here mean the sum of 
the squared differences between the A and the B image 
magnitudes (see below). 
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Fig. 5. Results from aperture photometry of QSO 0957+561A,B. There are 422 binned data points for each quasar image. The 
B data is shifted by —0.15 mag. Error bars are la limits. 



A; 



The data model is 



q(U) + e A (U) , i = 1, 



,N A 



Bj = q(tj - r) + %) + e B (tj) , j = 1, 



where q(t) denotes the inherent quasar variability, which 
should be the same in the two images. l(t) takes care of the 
unknown amplification ratio between A and B, as well as 
any long time scale microlensing. The observational errors 
are tA(ti) and esfe)- 

The combined light curve (denoted in the formulae as 
C) is constructed by taking the A values as they are and 
"correcting" the B data by l(t) and shifting them by r: 



A, 



if tk = U 



Bj - l(tj) if tk =tj +t 



(3) 



where k = 1, . . . , N and N = Na + Nb ■ The dispersion of 
the combined light curve (abbreviated as CLC in the text) 
is now computed for a range of r-values. The resulting 
dispersion spectra 



D\r) 



min D 2 (T,l(t)) 
i(t) 



(4) 



can subsequently be inspected with regard to minima. The 
time shifts of the most significant minima arc candidates 
for the true time delay. 

The accuracy of the observations is taken into account 
by using the statistical weights, Wi — 1/€a(U) and Wj 



l/esitj). The squared difference between two data points 
in the CLC (see estimates below) must be multiplied with 
the combined statistical weights Wk — Wij = w ' + w- ■ 

With Z(i)=constant, the A and B curves are considered 
to be unaffected by microlensing variability, differing only 
by the unknown ratio of the amplification factors of gravi- 
tational (macro-) lensing. We shall, however, also compute 
spectra where we account for slowly varying microlensing 
effects in one or both of the light curves. In these cases 
we set l(t) equal to a polynomial (typically of order two 
to eight). The B data is thus "modified" by the perturb- 
ing polynomial, into Bj + l{tj), and the coefficients of the 
polynomial are determined in such a way as to minimize 
the dispersion between A and B data. 

In this analysis we use two different methods to esti- 
mate dispersions. The simplest one is 



J2k=i SkWkGk(Ck+i — Ckf 
2 J2k=i SkWkGk 



= min 

i(t) 



(5) 



where Wk is the statistical weights, and Gk — 1 only if 
Cfe+i and Ck are from different data sets. (That is, one 
point from A and one from B. Otherwise, Gk = 0). Sk 
constrains the time gap between the AB or BA pairs; 



_ dk 



1 if 

if |tfc+i 



t k \ >S 



(G) 
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where 8 (measured in days) is the largest time separation 
allowed, also called the decorrelation length. 
The second statistic is 



D\ h = min 

' l(t) 



V JV-1 v jv Mh) w 



,.G 



,(7) 



where W n ^ m are statistical weights and G n . m controls that 
only AB or BA pairs are considered, just as in the previous 
estimate. In this scheme more pairs are included by not 
only considering strictly neighboring pairs; Snjn weights 
(C„ — C m ) 2 according to the corresponding temporal sep- 
aration, \t n — t m \. We may use a flat window (h = 1), 
linear (h = 2) or Lorcntzian (h = 3) down-weighting, see 
Pelt et al. (1996). Here we use linear down-weighting; 



5(2) 



K-t m \ 



if \tn 
if Hn 



< s 

> s 



(8) 



The D\ and D\ 2 estimates calculate the dispersion 
between A and B points, with an upper limit on the 
corresponding time separation. D\ includes only con- 
secutive AB (or BA) points, and does not involve any 
smoothing. The D\ 2 estimate has the advantage of in- 
cluding much more pairs, and thus suppressing noise in 
the dispersion spectra. However, one should be careful 
not to over-smooth the spectra by using large decorre- 
lation lengths. We shall employ the D\ 2 estimate most of 
the time. Different values of the decorrelation length, 5, 
will be tested, as well as various models (constant versus 
higher-order polynomials) to account for the additional 
(microlens-induced) variability. 

ISDA contains a simple bootstrap procedure for cal- 
culating the error bars for the minima in the dispersion 
spectra. The CLC is smoothed and the corresponding 
residuals (data points minus "smoothed curve values") 
are re-shuffled 1000 times to create bootstrap samples. 
Smoothing is performed using a 7-point median filter, with 
an upper limit (typically a few days) on the time separa- 
tion between two successive data points. 

The time delay value from a particular dispersion esti- 
mate is taken to be the mean of the time delay distribution 
(an example is given in Fig. 11). The standard deviation 
of the distribution gives the estimated error. We quote the 
la errors. 



5.1.1. Complete data set 

From the complete data set of 422 data points for each of 
the two quasar images, we discarded six outliers. 

The number of pairs included in the dispersion esti- 
mates depends on 8. Fig. 6 displays selected window func- 
tions for Z?| and D\ 2 . The window function is the num- 
ber of nearby AB (or BA) pairs in the CLC as a func- 
tion of time shift. Obviously, larger 8 yields more pairs 
in the computation, but the overall shape of the window 
functions remains more or less the same (the curves get 
smoother as 8 increases). The sampling of the observa- 
tions may disfavor some time shifts, i.e. the number of 



2000 



1 500 -A 



500 




300 
Shift (days) 

Fig. 6. Window functions for the complete data set when em- 
ploying the estimates D\ 2 {8 — 5 days; upper curve, 8 = 2 
days; middle curve) and D% , 5 = 5 days (lower curve) . 



pairs of nearby points in the CLC can be very low for cer- 
tain shifts. Fortunately, there are no major depressions in 
the curves, so the statistical reliability of the dispersion 
values should not vary much for the different trial shifts 
(especially in the interesting range 400-440 days) . This is 
a reassuring and important fact. 

We shall plot the dispersion spectra for trial shifts in 
the range 380-480 days. However, before we present and 
discuss the behavior of the dispersion curves in this limited 
range, we show in Fig. 7 a plot of two spectra calculated 
using the D\ 2 estimate where the interval goes from 0- 
600 days. Over the entire range, the dispersion is smaller 




100 



200 300 400 

Shift (doys) 



500 



600 



Fig. 7. D| 2 dispersion spectra, 8 = 3. The thin curve cor- 
responds to l(t) = lo, while the thick curve is computed by 
also accounting for additional fluctuations (l(t) is a 5th order 
polynomial). 



for the estimate that includes the perturbing polynomial. 
A higher-order polynomial would account even more for 
differences in the two quasar signals, and thus decrease the 
general dispersion even further. One must be wary not to 
"over-correct" the B data, though. 

We first computed dispersion spectra using various 
decorrelation lengths, but without any corrections for mi- 
crolensing (l(t) = Iq)- Then, the calculations were re- 
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peated, but this time we included polynomials to model 
long time scale ML variability in the light curves. The 
results were not very sensitive to the degree (2nd-8th or- 
der) of the perturbing polynomial. To limit the number of 
variables, we fixed the degree to 5th order. Figs. 8 and 9 
display a selection of spectra (5 = 1, 3, 5, 7 days) derived 
from the two methods. 
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Fig. 8. D\ 2 dispersion spectra with 8 = 1, 3, 5, 7 days. No 
corrections for additional, possibly microlens-induced variabil- 
ity. 
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Fig. 9. D| >2 dispersion spectra with 5 = 1, 3, 5, 7 days. A 5th 
order polynomial was used for modeling of additional variabil- 
ity in the light curves. 

The curves are smoother for larger 6. Note that 6=1 
hardly involves any smoothing. For a given S, the two 
methods (without/ with correction for additional variabil- 
ity) yield very similar results. In both cases, the position 
of the minimum increases slightly for increasing decorre- 
lation length, from 424 to 425.5 days. Secondary minima 
are almost always found around 410 days, but they are 
moderated as the smoothing increases. 

We also computed spectra for 5 = 1 — 20 with incre- 
ments of one, and noted the position of the corresponding 
minima. Fig. 10 plots shift values corresponding to min- 
imum dispersion as a function of the 6 parameter. Here 
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Fig. 10. Dependence of the 5 parameter for the statis- 
tic. Asterisks: no correction for additional variability. Open 
squares: 5th order polynomial. 

we see more clearly that the minima are shifted towards 
higher values as pairs with larger time separations are in- 
cluded in the estimates. However, for 5 > 14—15, the trend 
is reversed, but we are probably smoothing too much al- 
ready Because the sampling of the observations is rather 
good (85% of the time separations are less than five days), 
we can afford to use small decorrelation lengths. This is re- 
assuring, because it reduces the danger of bias from pairs 
with large time separations. A striking feature is also seen: 
whether we account for slowly varying microlensing effects 
or not, the minima are all at 424 days for 5 = 1, 2, 3 and 
4. 

The same procedure was employed for the D\ statis- 
tic, which only includes consecutive A and B points whose 
time separation in the CLC is less than 5. Hence, no 
smoothing is performed. These spectra were similar to the 
curves in Figs. 8 and 9 corresponding to 6 = 1. The results 
confirmed the trends and features which were highlighted 
above. Since the number of pairs included in this estimate 
is lower than the "smoothing" estimates, it is not as reli- 
able, statistically speaking. The spectra did exhibit more 
noise, but consistently produced global minima between 
423.5 and 424.5 days when S was < 4 days. As before, 
larger decorrelation lengths yielded higher time shift val- 
ues. 

To summarize: the different dispersion estimates (D|, 
D\ 2 , both with and without microlensing correction) all 
produce minima which are concentrated around 424 days 
for 5 < 4. It seems that larger decorrelation lengths intro- 
duce bias. Hence, we performed bootstrap runs only for 
the limited range of <5-values. As noted earlier, the results 
were not significantly affected by the nature of the poly- 
nomial used to model additional variability. In particular, 
for the preferred range of S (i.e. 1-4 days), the minima in 
the dispersion spectra were all at 424 days for degrees of 
order two to eight. We do not want to use very high-order 
polynomials, as this could suppress some of the intrinsic 
quasar fluctuations. Hence, it seems justified to fix the 
degree to five. 

Table 3 lists the results of the bootstrap procedure. In 
Fig. 11 we show an example of the distribution of time de- 
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Table 3. Time delay results for the D i 2 dispersion estimate 
applied to the complete aperture photometry data set. A 5th 
order polynomial was used to model additional variability in 
the quasar light curves. As noted in the text, other polynomials 
gave very similar results. Estimated errors are la limits. 



Statistic 


8 


Time delay (days) 


Dl, 2 


1 


424.8 ± 1.6 


D% 2 


2 


424.3 ± 0.6 


Dl, 2 


3 


424.8 ± 1.3 


Dl, 2 


4 


424.7 ± 0.7 



lays from one of the bootstrap runs. Here, the mode and 
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Fig. 11. Bootstrap results. The mean time delay from 1000 
bootstrap runs is indicated by the large tick mark. Statistic: 
-C*4,2! 8 = A, l(t) is a 5th order polynomial. 



the median were both 424.5 days, while the mean time 
delay was 424.7 days. The shapes of the other distribu- 
tions were similar, and they all had a small skew. We also 
note that the D\ estimate gave similar results with the 
bootstrap procedure. The estimated time delay using dif- 
ferent setups of the Dispersion estimation technique agree 
well. We thus take the most probable time delay to be the 
average of the numbers in the table, i.e. 424.7 days, with 
a mean estimated error of 1.1 days. However, we do not 
claim that the true error is as low as this. 

The magnitude difference, Attiab, between the A and 
B components in the 1992-1997 time span was ~ 0.076 
mag (we were not able to compute error bars for this pa- 
rameter). Hence, B was somewhat brighter than A in the 
time span which our data covers. This is commonly ex- 
plained by microlensing in the B component, see e.g. Pelt 
et al. (1998). 

Finally, it is worth noting that all spectra show a local 
maximum at ^418 days. So with this particular data set 
and the Dispersion estimation method, we can say that a 
time delay of roughly 418 (± 2) days seems rather unlikely. 

5.1.2. Truncated data sets 

We now estimate the time delay using only selected data 
segments. We divide each light curve into four parts, cor- 



responding to the seasons with J.D.— 2448000 roughly in 
the ranges 850-1200 (period 1), 1200-1600 (period 2), 
1600-1900 (period 3) and 1900-2300 (period 4), see Fig. 5. 
Assuming a time delay of around 425 days, it is clear that 
there is (fortunately) a large overlap between periods 1, 
2 and 3 of A and periods 2, 3 and 4 of B, respectively. 
We thus have the possibility of estimating the time delay 
between the quasar images from three different data sub- 
sets (let us call these 51, 52, 53). The motivation behind 
this is to see whether the truncated data sets all produce 
a minimum in the dispersion spectra around 425 days. We 
shall not perform an exhaustive analysis, though. 

Because the number of pairs included in the calcula- 
tions is much lower for these subsets, the statistical reli- 
ability is not as good as in the case where we used the 
complete data set. We employ the D| 2 estimate, 5 in the 
range 2-5 days, and compute spectra for trial shifts in 
the range 400-450 days. The effect of correcting for any 
non-intrinsic quasar fluctuations is also tested. For this we 
use (only) a 3rd order polynomial (the curves overlap for 
about 200 days in all three cases, and we do not allow for 
any extrinsic high-frequency components within this time 
span) . 

Fig. 12 displays D\ 2 dispersion spectra computed for 
the subsets 51, 52 and 53, assuming a constant magni- 
fication ratio, l(t) — Z . Allowing for a time-dependent 
magnification ratio (to account for possible microlcnsing- 
effects) did not significantly change the overall shape of 
the curves. The minima are sometimes split in two for 
short decorrelation lengths, so we use 8 — 4 which reduce 
the "noise". 

For the first data subset, 51, the deepest minimum 
is at 430 days. Secondary minima are seen around 424 
and 415 days. We checked the window function too, and 
it contained a prominent minimum for the 430 day shift. 
It might be that the 430 day candidate is caused by the 
combination of irregular sampling and an "unfortunate" 
time shift. The second subset, 52, yields two minima in 
the dispersion spectrum, 412 days and 425 days, the lat- 
ter being the deepest. Here, the window function had no 
unfavorable time shift. The third plot shows the results 
using the 53 data, and here the deepest minimum is posi- 
tioned around 425 days. Secondary minima are found for 
time shifts of 408 and 434 days. The window function had 
minima at 408 and 425 days, corresponding exactly to two 
of the observed local minima in the dispersion spectrum. 

Bootstrap runs were performed to get an estimate of 
the uncertainty. From the distribution of 1000 time delays, 
the results (mean and standard deviation) were as follows: 
430.6 ± 2.7 days (51), 424.9 ± 3.0 days (52), 426.1 ± 2.3 
(53). We do not attempt to judge the reliability of the 
different time delay candidates. It is interesting to see, 
however, that the local minima may be found in usually 
one of three regions, i.e. around 410, 425 and 430 days. 
Also, a local peak in the interval 416-420 days is found 
in all spectra, thus supporting the statement made in the 
previous section that the correct time delay is less likely 
to lie in this range. 
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Fig. 12. Z)| 2 dispersion spectra. 8 = 
4. No corrections for microlens-induced 
variability. The shapes of the disper- 
sion spectra from the different data sets 
are very different. The position of some 
minima correspond to minima in the 
window function, i.e. the number of 
pairs in the dispersion estimates reaches 
a minimum for these time shift values. 



For the complete data set, we found in the previous 
section that Atoab ~ 0.076 mag. However, a brief look at 
the combined light curve (B shifted in time by —425 days 
and in magnitudes by 0.076 mag) indicated that a single 
magnitude difference did not optimally align the A and 
B data. We thus computed Atoab for each of the three 
data subsets in order to investigate this further. We got 
i=a 0.089 mag, 0.086 mag and w 0.050 mag for subsets 
ST, 52 and 53, respectively. The fact that the magnitude 
difference seems to be time-dependent will be discussed in 
Sect. 6. 

5.1.3. Time delay from RES's data 

Now that we have performed a time delay determination 
using our new photometric results, it might be interesting 
to see whether the "old" reductions by Schild and collab- 
orators 5 give similar results. We used 537 data points for 
each quasar image which covered the same period as our 
data. An extensive analysis is beyond the scope of this 
paper. We shall thus only comment on the main results 
from the D\ 2 dispersion estimate. 

The results revealed a few interesting general trends: 

— With short decorrelation lengths (5 < 5 days) and no 
correction for microlensing, the dispersion spectra had 
minima at 435 days. Prominent, but marginally higher 
minima were seen at ~412 days. By introducing poly- 
nomials, l(t), of varying degrees to model (hypotheti- 
cal) microlensing, the global minimum was shifted to 
around 411 days for the same range of S. 

— For larger decorrelation lengths, the minima were 
found in the region 411-415 days, irrespective of the 
exact nature of l(t) (constant versus higher-order poly- 
nomials). Increasing S from 6 to 12 days consistently 
produced minima at larger time shifts, going from 412 
up to 415 days. 

We present in Fig. 13 dispersion spectra D\ 2 with 5=1, 
3, 5 and 7 days. 

5 See data table at http://cfa-www.harvard.edu/~rschild/. 
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Fig. 13. Dispersion spectra D\ 2 = 1, 3, 5, 7 days) with 5th 
order polynomial, using RES's data. 

The second point above describes a general trend seen 
in all the different analyzes : the position of the dispersion 
minimum increases with increasing decorrelation length. 

Here we find that the position of the global minimum 
was 411 days for 5 < 4 days, and it did not depend on the 
order of the perturbing polynomial (2nd to 8th order). 
There are certainly indications of additional variability in 
the observational data, hence it seems justified to intro- 
duce the l(t) polynomial into the dispersion minimization 
process. The results of the bootstrap procedure are listed 
in Table 4. The mean time delay is 411.7 days, with a 
mean estimated error of 1.9 days. 

5.2. x 2 minimization 

The method of Burud et al. (2001) is based on \ 2 min- 
imization between the data and a numerical model light 
curve. Microlensing in one or both light curves may be cor- 
rected for. We shall in the following carry out a brief, non- 
exhaustive time delay analysis with this method. Because 
the procedure is explained in detail in the paper by Burud 
et al., we only summarize the main features. 
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Table 4. Time delay results for the D\ 2 dispersion estimate 
applied to the RES photometry data. A 5th order polynomial 
was used to model additional variability in the quasar light 
curves. Estimated errors are la limits. 



Statistic 


8 


Time delay (days) 


Di, 2 


1 


411.4 ± 2.3 


D\ 2 


2 


411.8 ± 1.8 


Di, 2 


3 


411.8 ± 1.8 


D\ 2 


4 


411.8 ± 1.7 



The underlying idea is that, in the absence of ML, one 
can model the two quasar light curves with one model 
curve, g(t), together with two parameters, r and Am, 
describing the time shift and magnitude offset between 
images A and B. An arbitrary model curve with equally 
spaced sampling points is x 2 minimized to the two orig- 
inal light curves. The minimization is done only for the 
observed data points, so only the model curve is interpo- 
lated and not the data. Because the data is irregularly 
sampled and contain noise, a smoothing scheme is neces- 
sary. Here, the model light curve, g(t), is smoothed on a 
time scale, 7\, corresponding to the typical sampling in- 
terval of the data. The smoothing term is multiplied by 
a Lagrange parameter, A, which can be chosen so that 
the model curve matches the data correctly in a statisti- 
cal sense for adopted Gaussian statistics (this parameter 
has no physical meaning). In addition, each data point is 
given a weight which depends on the relative distances to 
all other points in the curve. Down-weighting is performed 
using a Gaussian with FWHM = 2V21n2T 2 , where the 
user may choose the T 2 parameter. The weights are nor- 
malized, so that the maximum value of Wi is 1. This would 
be the case if only one point is within the time interval 
defined by T 2 . According to the authors, a sensible choice 
of T 2 would be the approximate time scale of the intrinsic 
quasar fluctuations. 

The method was only applied to the complete light 
curves. We did not attempt to model higher-order ML 
fluctuations in the light curves, as this is quite an elaborate 
process. A wide range of parameter setups was tested, and 
the results proved to be remarkably stable. We present 
only the main results. 

The x 2 -values as a function of time shift (in the range 
400-450 days) are plotted in Fig. 14. The parameters were 
as follows: Ti = 4 days, T 2 = 20 days and A = 4000. We 
recognize some of the features from the analysis with the 
Dispersion estimation technique: The minimum x 2 -value 
occurs for a time shift of 425 days. A secondary minimum 
is found around 411 days, but the x 2 is n °t as low as the 
tiny, local minimum at 431-432 days. The overall shape 
of the x 2 distribution remains the same even for large 
variations in the parameters. The lowest x 2 -value is always 
in the range 424-426 days. 

The results form Monte Carlo simulations yielded a 
time delay of 425.1 ± 1.3 days. Also with this method 
we find that a time delay of roughly 415-420 days seems 
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Fig. 14. x as a function of time shift. 



less likely; the x 2_ curve typically has a maximum in this 
range. 



6. Microlensing 

We will now briefly investigate the microlensing residual 
in the quasar light curves. The standard procedure is to 
shift the light curves in time to correct for the different 
light travel times, and then subtract them from each other. 
The last step is not trivial, as the A and B data points are 
irregularly sampled. This means that when we shift the 
B data in time by — r, the A and B points will generally 
not overlap. The B data point to be subtracted from a 
particular A point might be several days away. We have 
addressed this issue in a simple way. 

After having shifted the B curve by —425 days, we 
check for each data point of image A whether there is 
a point from the B curve within a certain gap limit of 
the current A point. If this is the case, then B is sub- 
tracted from A, and the result is stored in a "residuals 
array" with an averaged time argument. The procedure 
only makes use of the original, raw data points, and there 
is one free parameter, namely the gap limit. The exact 
value of this parameter depends mostly on the spacing of 
the data, but also on the (assumed) ML time scale. With 
good sampling the gap limit can be set quite low (a few 
days) and thus, at least in principle, enable investigation 
of rapid ML. Lowering the gap limit will obviously de- 
crease the number of points in the residuals array. On the 
other hand, the A — B differences are then calculated from 
AB pairs with smaller time separations, which is a good 
thing if one wants to probe short time scale fluctuations 
in the residuals. 

With a time delay of ~425 days, our A and B data 
overlap for about 3.5 years. The ML investigation will 
consequently only cover this time span. The large tem- 
poral gaps in the residuals data (a result of the lack of 
observations in the months where 0957+561 was below 
the horizon) also precludes a continuous "signal" for the 
whole period. 
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6.1. Trends on long time scales 

In Fig. 15 we display the A — B residuals (Attiab), com- 
puted by following the above procedure. We adopted a 
time delay of 425 days, and the gap limit was set to 2.5 
days. There are three seasons which contain an adequate 
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Fig. 15. Time-shifted B curve values subtracted from A curve 
values. Dates on the abscissa relate to the (unshifted) A 
light curve. The A — B differences, AmAB, are clearly time- 
dependent. Three examples of spline approximations to the ML 
residuals are also shown. Red, green and blue lines correspond 
to 6, 9, and 11 nodes, respectively. 

number of points. The results look rather noisy, but there 
are some significant features. In particular, the third sea- 
son clearly has a different magnification ratio between A 
and B, compared to the first two. Moreover, it also varies 
within the particular season. Variability on shorter time 
scales may also be seen at certain periods. 

The amplitude of the variations in the first two seasons 
is about 0.05 mag. From the plot we can see that the av- 
erage magnitude offset for the two first seasons are ^0.09 
mag, while for the third period the number is roughly 
0.05 mag. This is in good agreement with the values we 
obtained in Sect. 5.1.2. 

Here, we only want to get an idea of the general trends 
in the ML residuals. We thus tried fitting standard cubic 
splines with different number of nodes into the residuals, 
see Fig. 15. The (red) curve with 6 nodes "detects" only 
changes on long time scales. It is thus a rather conservative 
guess as to how microlenses in the macro-lensing galaxy af- 
fected the light from the quasar images. The general trends 
are similar to the results of Pelt et al. (1998) - compare 
with the three last "seasons" in their Fig. 9. Pelt et al. 
do not have the points around J.D.-2448000=2100 which 
we do (see Fig. 15). Although very sparse, these data in- 
dicate that the curve does not continue to fall off. The 
other two splines probe finer details in the ML residuals, 
and are more optimistic approximations. Some of the vari- 
ability, notably in the gaps, is highly questionable. Still, 
the ~0.05 mag drop in the difference data from the sec- 
ond to the third season (on the order of 300 days) seems 
significant. We conclude that microlensing variability of 



approximately 5% amplitude on time scales of less than a 
year has been significantly observed. 

6.2. Short time scale variation 

On short time scales (the order of weeks) the residuals 
are rather noisy. Fig. 16 shows the three first seasons of 
the ML residuals in greater detail. For two of the sea- 
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Fig. 16. The ML residuals from the three first seasons of 
Fig. 15. For two periods we have fitted a cubic spline (upper 
panel; 9 nodes, lower panel; 4 nodes). 



sons we also include spline approximations. As can be seen 
from the middle panel, the data here shows no significant 
trends, thus no interpolations are attempted. (Compare 
the scatter here to the the mean formal errors in the aper- 
ture photometry of the quasar images, i.e. 17 mmag). It 
is not trivial to extract information on the true microlens- 
induccd fluctuations from data such as these, and the 
interpolated curves are mostly meant to guide the eye. 
However, on the first plot we can discern a steep negative 
slope in the residuals around J.D. -2448000 = 900 followed 
by a positive slope some 50 days later. Optimistically, we 
can explain this as an "event" lasting on the order of 70 
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days, but nothing certain can be said about its amplitude. 
It could potentially be a strong event, because of the steep 
gradients. The third period (lower panel) indicates more 
clearly a U-shaped feature. The amplitude and time scale 
is hard to assess, because we do not have any points in 
the "wings". But is seems that the amplitude is at least 
0.05 mag, and the time scale is on the order of 200 days. 

7. Summary and discussion 

We have presented a re-reduction of the CCD image 
frames critical to the discussion about time delay and 
microlensing in the 0957+561 gravitational lens system. 
Improved computational techniques allow better subtrac- 
tion of the effects of the lens galaxy, and correction for the 
aperture crosstalk that arises in aperture photometry of 
the somewhat overlapping quasar images. 

Analysis of the re-reduced photometry for time delay, 
principally using several variants of the Dispersion tech- 
nique, gives consistent values around 425 days. The aver- 
age result from the Dispersion method and the \ 2 mini- 
mization method is 424.9 days, with an estimated mean 
error of 1.2 days. However, we do not claim that the true 
error is as small as this. We also note that time delays 
of roughly 416 to 420 days were never seen in this in- 
vestigation and are thus less favored by us. This is not in 
agreement with e.g. Kundic et al. (1997), Pelt et al. (1998) 
and Colley & Schild (2000). 

Analysis of principally the same image frames with 
fundamentally different reduction and time delay estima- 
tion techniques had previously given 404 days (Schild & 
Thomson 1997, Direct Autocorrelation) and 416.3 days 
(Pelt et al. 1998, Dispersion estimation procedure), but 
re-analysis by Oscoz et al. (2001) of the same bright- 
ness record gave estimates near 422.6 days. Other smaller 
data sets for approximately the same observational epochs 
gave 417 days (Kundic et al. 1997, PRH method, Linear 
Interpolation) and 425 days (Serra-Ricart et al. 1999, S 2 
method). 

In all cases but the first, the quoted errors (typically 
a day or two) are much smaller than the discrepancies 
between different data sets or between estimates for the 
same brightness record. Critical to the discussion is the 
fact illustrated in Fig. 7 that the FWHM of the disper- 
sion curve has a value of approximately 100 days, and even 
the local minima, as seen in Figs. 8, 9, 12 and 14, have 
a FWHM of 10-20 days in spite of the daily data sam- 
pling and the available several hundred data points for 
any test lag (Fig. 6). It is by now evident that something 
fundamental limits our ability to estimate time delay to 
the expected limits imposed by the data sampling and the 
observational errors. 

The physical origin of this discrepancy has been at- 
tributed by Colley et al. (2003) to the fact that the 
quasar's luminous structure is time-resolved and mi- 
crolensed. This combination of microlensing and a time- 
resolved source might produce multiple time delays whose 
pattern changes from year to year. Some evidence for this 



may be seen in Fig. 12, where the relative importance of 
persistent lags near 410, 425, and 430 days seems to have 
changed during the observational period. We note, how- 
ever, that for some subsets these time lags coincide with 
minima in the window functions. 

This puts a new perspective to the understanding of 
the role of microlensing for a quasar source. Previous 
discussions of microlensing (see Schmidt & Wambsganss 
1998, and references contained therein) have focused upon 
the role of small accretion discs crossing the network 
of caustics in the magnification diagram produced by 
MACHOs in the lens galaxy. If real, microlensing fluc- 
tuations on time scales on the order of 70 days (Sect. 6.2) 
may signal the presence of MACHOs with masses possi- 
bly down to planetary masses. The small, 5% amplitude 
of this short time scale microlensing signal is consistent 
with previous conclusions that the luminous source may 
be quite large relative to the Einstein Rings (Refsdal & 
Stabell 1991, 1993, 1997; Refsdal et al. 2000). Because 
this is near to the noise level, extremely careful data ac- 
quisition and analysis is called for in determining the time 
delay and microlensing. 

We hope in future papers to extend the time delay 
and microlensing analysis, using an even larger data set. 
A longer observational base line and maybe more statisti- 
cal techniques could shed new light on the time delay issue. 
We also hope that our new reduction scheme (both aper- 
ture and PSF photometry, see Ovaldsen 2002) includes 
some new features which could be of interest to other re- 
searchers. 
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